---
title: "Government Religious Discrimination, Support of Religion, and Muslim
Minority-Related Violence in Western Democracies"
output: html_documentR.version
date: "2023-11-24"
---

*Note: The analyses were conducted using R version 4.3.1. All necessary packages are indicated where applicable. The corresponding labels for the tables and figures in both the paper and the appendix are provided at the start of the code snippets.


*Descriptives*


```{r}
### Descriptive statistics (Table. 1) 

load("final_dataset.rdata")

Descriptive_main <- subset_data <- final_data[, c("WZMIN2MAJX","MMXX","ln_population","ln_GDP", "YouthRatio.y","LXX","soc_vioc", "MINPCTX", "gini_disp",  "logland",  "popvote_im" ,"church", "relimp", "majority", "bur_qua")] 

Descriptive_main <- as.data.frame(Descriptive_main)

library(stargazer)
stargazer(Descriptive_main, title = "Descriptive Summary Table", digits=1, out="table3.txt",type = "text",
          covariate.labels=c("Violence by Muslims", "DGRD against Muslims","Population","GDP Per Capita", "Youth Bulge",
                              "Support for Religion", "Violence by Christian Groups", "Muslim Minority Percent", "Economic Inequality", "Land Area(logged)",  "Populist Party Support", "Religious Attendance", "Religious Importance", "Majority Group Percent",  "Bureaucratic Quality"))


stargazer(Descriptive_main)

```





```{r}
descriptive_M <-  na.omit(unique(subset(final_data, select=c(WZMIN2MAJX, soc_vioc, cyear, year,country, MMXX, LXX,MMXX_d, lag_MMXXd))))
  
### Within-Between variation of the chief dependent and independent variables reported in the paper (P.18-19)


# GRD

library(dplyr)
  
descriptive_M <-descriptive_M %>%
  group_by(country) %>%
  mutate(N_changesMX = sum(diff(MMXX) != 0)) 

descriptive_M <-   descriptive_M %>%
  group_by(country) %>%
  mutate(country_meanMX = mean(MMXX))

Within_MX <- descriptive_M %>% group_by(country) %>% summarise(within_changes = sum(unique(N_changesMX)))
sum(Within_MX$within_changes)
length(unique(descriptive_M$country_meanMX))


# Minority violence

descriptive_M <- descriptive_M %>%
  group_by(country) %>%
  mutate(N_changesminvio = sum(diff(WZMIN2MAJX) != 0))

# Mean
descriptive_M <-   descriptive_M %>%
  group_by(country) %>%
  mutate(country_meanminvio = mean(WZMIN2MAJX))


Within_minvio <- descriptive_M %>% group_by(country) %>% summarise(within_changes = sum(unique(N_changesminvio)))
sum(Within_minvio$within_changes)
length(unique(descriptive_M$country_meanminvio))


# Majority violence

descriptive_M <- descriptive_M %>%
  group_by(country) %>%
  mutate(N_changessocvio = sum(diff(soc_vioc) != 0))

# Mean

descriptive_M <-   descriptive_M %>%
  group_by(country) %>%
  mutate(country_meansocvio = mean(soc_vioc))

Within_socvio <- descriptive_M %>% group_by(country) %>% summarise(within_changes = sum(unique(N_changessocvio)))

sum(Within_socvio$within_changes)

length(unique(descriptive_M$country_meansocvio))

# GSR

descriptive_M <- descriptive_M %>%
  group_by(country) %>%
  mutate(N_changesLXX = sum(diff(LXX) != 0))

# Mean
descriptive_M <-   descriptive_M %>%
  group_by(country) %>%
  mutate(country_meanLXX= mean(LXX))

Within_LXX <- descriptive_M %>% group_by(country) %>% summarise(within_changes = sum(unique(N_changesLXX)))
sum(Within_LXX$within_changes)

length(unique(descriptive_M$country_meanLXX))

```
*Descriptive Graphs*

```{r}
### Mean graph of all countries for each year (Figure.1)

# Calculate the mean of each variable for all countries
means_1 <- descriptive_M %>%
  group_by(year) %>%
  summarise(mean_min = mean(WZMIN2MAJX))

means_2 <- descriptive_M%>%
  group_by(year) %>%
  summarise(mean_soc = mean(soc_vioc))

means_3 <- descriptive_M%>%
  group_by(year) %>%
  summarise(mean_mx = mean(MMXX))

means_4 <- descriptive_M%>%
  group_by(year) %>%
  summarise(mean_lxx = mean(LXX))

means <-  merge(means_1, means_2)

means <- merge(means, means_3)

means <- merge(means, means_4)



# Majority violence 

library(ggplot2)
d1 <- ggplot(means, aes(x = year, y = mean_soc, color = "Mean Soc")) +
  geom_line(show.legend = FALSE) +
  scale_color_manual(values = "sienna") +
  labs(title = "",
       x = "",
       y = "Majority Violence (Mean)") +
  theme_minimal()




# Minority violence 

d2 <- ggplot(means, aes(x = year, y = mean_min, color = "Mean min")) +
  geom_line(show.legend = FALSE) +
  scale_color_manual(values = "azure4") +
  labs(title = "",
       x = "",
       y = "Muslim Violence (Mean)") +
  theme_minimal()



#  GRD

d3 <- ggplot(means, aes(x = year, y = mean_mx, color = "Mean mxx")) +
  geom_line(show.legend = FALSE) +
  scale_color_manual(values = "lightblue") +
  labs(title = "",
       x = "Year",
       y = "DGRD against Muslims(Mean)") +
  theme_minimal()


# GSR

d4 <- ggplot(means, aes(x = year, y = mean_lxx, color = "Mean lxx")) +
  geom_line(show.legend = FALSE) +
  scale_color_manual(values = "lightsalmon") +
  labs(title = "",
       x = "Year",
       y = "Support for Religion (Mean)") +
   ylim(8, 10) + 
  theme_minimal() 


# save descriptives 

library(gridExtra)

descriptives <- grid.arrange(d1, d2, d3, d4, ncol=2, nrow=2)



### Within-country variation (Figures A1-A4 in the appendix)


## Minority Violence (Figure.A1)

#Select countries where within variation exists

selected_countries <- c("Australia", "Belgium", "Canada", "Greece", "France","Denmark","Netherlands", "Norway","Sweden","Switzerland", "United Kingdom", "United States of America")

# Filter the data to include only the selected countries
descriptive_filt <- descriptive_M[descriptive_M$country %in% selected_countries,]

within_bycountry <- ggplot(descriptive_filt, aes(x = year, y = WZMIN2MAJX, group = country, color = country)) +
  geom_line(show.legend = FALSE) +
  facet_wrap(~ country,  scales = "free_y") +
  labs(title = "",
       x = "Year",
       y = "Muslim Minority Violence") +
  theme_minimal()


## Majority Violence (Figure.A2)

exclude_countries <- c("Portugal", "Finland")
descriptive_maj <- descriptive_M[!(descriptive_M$country %in% exclude_countries),]


within_maj <- ggplot(descriptive_maj, aes(x = year, y = soc_vioc, group = country, color = country)) +
  geom_line(show.legend = FALSE) +
  facet_wrap(~ country,  scales = "free_y") +
  labs(title = "",
       x = "Year",
       y = "Majority Group Violence") +
  theme_minimal()


## GRD (Figure.A3)

exclude_countries <- c("Portugal", "Ireland","Liechtenstein", "Luxembourg", "Canada","Australia")
descriptive_mx <- descriptive_M[!(descriptive_M$country %in% exclude_countries),]



within_dgrd <- ggplot(descriptive_mx, aes(x = year, y = MMXX, group = country, color = country)) +
  geom_line(show.legend = FALSE) +
  facet_wrap(~ country,  scales = "free_y") +
  labs(title = "",
       x = "Year",
       y = "GRD") +
  theme_minimal()


## GSR (Figure.A4)

exclude_countries_lxx <- c("Austria" ,"Liechtenstein", "New Zealand", "Netherlands", "France", "Belgium", "Denmark")
descriptive_lxx <- descriptive_M[!(descriptive_M$country %in% exclude_countries_lxx),]


within_lxx <- ggplot(descriptive_lxx, aes(x = year, y = LXX, group = country, color = country)) +
  geom_line(show.legend = FALSE) +
  facet_wrap(~ country,  scales = "free_y") +
  labs(title = "",
       x = "Year",
       y = "Government Support of Religion") +
  theme_minimal()


#### Means of variables for each country (Figure.A5)


# Minvio
library(countrycode)
descriptive_M$country_abr <- countrycode(descriptive_M$country, origin = "country.name", destination = "iso2c")

#exclude coutnries with 0 values

descriptive_minmean <- descriptive_M[descriptive_M$country_meanminvio!= 0,]


mean_min <-  ggplot(descriptive_minmean, aes(x = reorder(country_abr,country_meanminvio), y = country_meanminvio)) +
  geom_bar(stat = "identity", fill = "azure4",position = "dodge") +
  #geom_text(aes(label = country_meanminvio), vjust = -0.5, position = position_dodge(width = 0.9)) +
  labs(title = "",
       x = "",
       y = "Muslim Minority Violence") + 
  ylim(0,2.5) + 
  theme_minimal()
  
  
# Majvio
  
  descriptive_majmean <- descriptive_M[descriptive_M$country_meansocvio!= 0,]  
  
mean_maj <-  ggplot(descriptive_majmean, aes(x =  reorder(country_abr,country_meansocvio), y = country_meansocvio)) +
  geom_bar(stat = "identity", fill = "sienna",position = "dodge") +
 # geom_text(aes(label = country_meansocvio), vjust = -0.5, position = position_dodge(width = 0.9)) +
  labs(title = "",
       x = "",
       y = "Majority Group Violence") + 
  ylim(0,10) + 
  theme_minimal()   
  

#  DGRD
  
descriptive_meanmx <- descriptive_M[descriptive_M$country_meanMX != 0,]    
  
mean_dgrd <-  ggplot(descriptive_meanmx, aes(x =  reorder(country_abr,country_meanMX), y = country_meanMX)) +
  geom_bar(stat = "identity", fill = "lightblue",position = "dodge") +
  #geom_text(aes(label = country_meanMX), vjust = -0.5, position = position_dodge(width = 0.9)) +
  labs(title = "",
       x = "Country",
       y = "GRD") + 
  theme_minimal()  
  
  
# GSR
 
mean_lxx <-   ggplot(descriptive_M, aes(x =  reorder(country_abr,country_meanLXX), y = country_meanLXX)) +
  geom_bar(stat = "identity", fill = "lightsalmon",position = "dodge") +
 # geom_text(aes(label = country_meanLXX), vjust = -0.5, position = position_dodge(width = 0.9)) +
  labs(title = "",
       x = "Country",
       y = "Government Support of Religion") + 
  theme_minimal() 
  
 
 # save all as one 

library(gridExtra)
  
mean_4 <- arrangeGrob(mean_min, mean_maj, mean_dgrd,mean_lxx,  ncol=2, nrow=2)


```


*REWB Models*

*Muslim violence (Table. 2)*
```{r}
# Null model (M0)
library(lme4)
library(lmtest)

Min0 <- lmer(WZMIN2MAJX~ 1 + (1 | country), final_data)
summary(Min0)


# Only controls (M1)


Min1 <-lmer(WZMIN2MAJX~  MINPCTX+ gini_disp +YouthRatio.y + logland +ln_population + ln_GDP + time +  (1 | country), final_data) 

summary(Min1)


# controls + GRD (M2)

Min2 <-lmer(WZMIN2MAJX~  MINPCTX+ gini_disp +YouthRatio.y + logland +ln_population + ln_GDP + time + lag_MMXXd+ + lag_MMXXmean_int + (1 | country), final_data)

summary(Min2)


# Full model(M3)

Min3 <-lmer(WZMIN2MAJX~  MINPCTX+ gini_disp +YouthRatio.y + logland +ln_population + ln_GDP + time + lag_MMXXd+ lag_MMXXmean_int +lag_lxxd_int + lag_lxxmean_int+ (1 | country), final_data)


summary(Min3)

# Compare min3 and min2 models

lrtest(Min2, Min3)


# Plot marginal effect (Figure 2)

library(sjPlot)

plot_model(
 Min3,
  type = "pred", 
 pred.type = "fe",
  terms = c("lag_MMXXd"), 
  title = "",
  ci.lvl = 0.95,
 auto.label = F,
  axis.textsize = 10,
  color = "lightblue",
  line.color = "black",  # Set line color
  line.size = 1.5, 
 size = 3) + 
theme_minimal() + ylim(-0.5, 1)  + labs(x = "GRD against Muslims", y = "Muslim Minotiry Violence")




```


*Majortiy Violence (Table. 3)*


```{r}
# Null model (M4)


Maj0 <- lmer(soc_vioc_int~ 1 + (1 | country), final_data)
summary(Maj0)


# Controls (M5)

Maj1 <-lmer(soc_vioc_int~ MINPCTX + YouthRatio.y+ logland+ ln_population + ln_GDP + popvote_im + church + relimp + majority + bur_qua +  time + (1 | country), final_data) 

summary(Maj1)


# controls + GRD (M6) 

Maj2 <-update(Maj1, . ~ .+lag_MMXXd + lag_MMXXmean_int ) 


summary(Maj2)


# Full Model (M7)

Maj3 <-update(Maj2, . ~ .+ lag_lxxd_int +lag_lxxmean_int)


summary(Maj3)


lrtest(Maj2,Maj3)


# Plot marginal effect (Figure.3)

 
plot_model(
 Maj3,
  type = "pred", 
 pred.type = "fe",
  terms = c("lag_MMXXd"), 
  title = "",
  ci.lvl = 0.95,
 auto.label = F,
  axis.textsize = 10,
  color = "coral",
  line.color = "black",  # Set line color
  line.size = 1.5, 
 size = 3) + 
theme_minimal() + ylim( 0 , 5) + labs(x = "GRD against Muslims", y = "Violence by Majority Groups")

```

*ROBUSTNESS Checks*

```{r}
# 1-  Extra control variables (Table A1-appendix)

Min4 <- update(Min3, . ~ . +  UnempTotILO.x + Religion + transnationalterrorism +educyears )

summary(Min4)

# majority

Maj4 <- update(Maj3, . ~ . +  UnempTotILO.x +Religion  +  transnationalterrorism +  educyears)


summary(Maj4)


# 2 - Sensitivity analyses using Sensemakr (Figure A11) 


library(sensemakr)

sensemakr_min3 <- lm (WZMIN2MAJX_int~ lag_MMXX+ lag(LXX) + MINPCTX+ gini_disp +YouthRatio.y + logland +ln_population  + ln_GDP  + UnempTotILO.x + Religion + transnationalterrorism +  factor(year) + factor(country), data =final_data )


#Benchmark: Minority Percent
 
sensitivity_min1 <- sensemakr(model=sensemakr_min3, 
          treatment=  "lag_MMXX", 
          benchmark_covariates = c("MINPCTX"),
          kd =1:2) 
 


plot(sensitivity_min1)


#Benchmark: econ_inequality

sensitivity_min <- sensemakr(model=sensemakr_min3, 
          treatment=  "lag_MMXX", 
          benchmark_covariates = "gini_disp",
          kd =3) 



plot(sensitivity_min)


# Majority

sensemakr_mmaj3 <- lm ( soc_vioc_int~ lag_MMXX+  lag(LXX) +ln_population + ln_GDP + popvote_im + church + relimp + majority +MINPCTX + bur_qua + YouthRatio.y + logland+ UnempTotILO.x + Religion + transnationalterrorism + factor(year) + factor(country), data =final_data )


# Benchmark: Minority Percent
  
sensitivity_maj1 <- sensemakr(model=sensemakr_mmaj3, 
          treatment=  "lag_MMXX", 
          benchmark_covariates = "MINPCTX",
          kd =1:2) 



plot(sensitivity_maj1)


# Benchmark: Majority Percent

sensitivity_maj2 <- sensemakr(model=sensemakr_mmaj3, 
          treatment=  "lag_MMXX", 
          benchmark_covariates = "majority",
          kd = c(3)) 
 

plot(sensitivity_maj2)



# 3 - Robust standart erros (Table A2)

# Muslim violence
library(clubSandwich)

rse_min3 <- vcovCR(Min3, cluster = final_data$country,  type = "CR0")

coef_table <- coef(summary(Min3))

coef_table <- as.data.frame(coef_table)

coef_table$robust_se <- sqrt(diag(rse_min3)) # Coefs and robust standard errors as table
print(coef_table)


# Majority violence  (Table A2 con't)

library(robustlmm)

Maj4_rob <-rlmer(soc_vioc_int~  lag_MMXXd  + lag_lxxd_int +lag_MMXXmean_int +lag_lxxmean_int + time + logland +ln_population + ln_GDP + popvote_im + church + relimp + majority +MINPCTX + bur_qua + YouthRatio.y + (1  | country), final_data) 

summary(Maj4_rob)


# 4- Granger Causality

gt_maj0 <- lmer(soc_vioc_int~  lag_socvioc + (1 | country), final_data)

gt_maj1 <- lmer(soc_vioc_int~  lag_socvioc + lag_MMXX + (1 | country), final_data)


lrtest(gt_maj0, gt_maj1)



gtr_maj0 <- lmer(MMXX~ lag_MMXX  + (1 | country), 
final_data)

gtr_maj1 <- lmer(MMXX~  lag_MMXX +lag_socvioc   + (1 | country), 
final_data)

lrtest(gtr_maj0, gtr_maj1)


# Granger test for  Minority


gt_min0 <- lmer(WZMIN2MAJX~  lag (WZMIN2MAJX) +  (1 | country), final_data)

gt_min1 <- lmer(WZMIN2MAJX~ lag(WZMIN2MAJX) + lag_MMXX+ (1 | country), final_data)


lrtest(gt_min0, gt_min1)

gtr_min0 <- lmer(MMXX ~  lag_MMXX +  (1 | country), final_data)

summary(gtr_maj1)

gtr_min1 <- lmer(MMXX ~  lag_MMXX +lag_WZMIN2MAJX +   (1 | country), final_data)

lrtest(gtr_min0, gtr_min1)


# 5 - # Muslim minority Model (M3) with GLM  (Table A3)

min3_glm <- glmer(minbin~  MINPCTX+ gini_disp +YouthRatio.y + logland +ln_population + ln_GDP  + lag_MMXXd+ lag_MMXXmean_int +lag_lxxd_int + lag_lxxmean_int+factor(time) + (1 |country), final_data, family=binomial, nAGQ=0)

summary(min3_glm)

# 6 - Models without outlier (Table A3)

## Minority 

# Exclude Outliers
library(HLMdiag)
res <- hlm_resid(Min3,include.ls = FALSE)
res2 <-  hlm_resid(Min3, level = "country", include.ls = FALSE)

residuals <- res$.resid

# Calculate z-scores for the residuals
z_scores <- abs(scale(residuals))

# Define a threshold for outliers
outliers <- z_scores > 2

outlier_observations <- which(outliers)

# Create a new dataset excluding the identified outliers
clean_outlier2<- final_data[-outlier_observations, ]

Min3_outlier2 <-lmer(WZMIN2MAJX~  MINPCTX+ gini_disp +YouthRatio.y + logland +ln_population + ln_GDP  + lag_MMXXd+ lag_MMXXmean_int +lag_lxxd_int + lag_lxxmean_int + time +  (1 | country), clean_outlier2)

summary(Min3_outlier2)

## Majority

res_maj <- hlm_resid(Maj3,include.ls = FALSE)
res2maj <-  hlm_resid(Maj3, level = "country", include.ls = FALSE)

residuals_maj <- res_maj$.resid
Q1 <- quantile(residuals_maj, 0.25)
Q3 <- quantile(residuals_maj, 0.75)
IQR <- Q3 - Q1

# Define lower and upper bounds
lower_bound <- Q1 - 1.5 * IQR
upper_bound <- Q3 + 1.5 * IQR

outliers_maj <- which(residuals_maj < lower_bound | residuals_maj > upper_bound)

# Create a new dataset excluding the identified outliers
clean_outlier2_maj<- final_data[-outliers_maj, ]

maj4_outlier2 <-lmer(soc_vioc_int~MINPCTX + YouthRatio.y+ logland+ ln_population + ln_GDP+lag_MMXXd+ lag_MMXXmean_int +lag_lxxd_int + lag_lxxmean_int +time+ popvote_im + church + relimp + majority + bur_qua + YouthRatio.y +  (1  | country) , clean_outlier2_maj) 

summary(maj4_outlier2)


# 7 within-between effect of control variables (Table. A4)

Min3_seperated <-lmer(WZMIN2MAJX~  lag_MMXXd  +  lag_MMXXmean_int + lag_lxxd_int +lag_lxxmean_int + minpercent_d + minpercent_m +  YouthRatio_d + YouthRatio_m + ln_GDP +  logland  + ln_population +time +gini_disp_d +gini_disp_m+    (1 | country), final_data) 

summary(Min3_seperated)


Maj3_seperated <-lmer(soc_vioc_int~  lag_MMXXd  +  lag_MMXXmean_int + lag_lxxd_int +lag_lxxmean_int + minpercent_d + minpercent_m  +  YouthRatio_d  +YouthRatio_m +ln_GDP  +  logland + ln_population + time+ majority_d + majority_m +  relimp_d + relimp_m + church_d + church_m +  bur_qua + popvote_im+  (1 | country), final_data)

summary (Maj3_seperated)

tab_model(Min3_seperated, Maj3_seperated)


# 8-Models with narrow definition of violence (Table. A6)


Min3vio_only <-lmer(minvio_only~  MINPCTX+ gini_disp +YouthRatio.y + logland +ln_population + ln_GDP + time + lag_MMXXd+ lag_MMXXmean_int +lag_lxxd_int + lag_lxxmean_int+ (1 | country), final_data)

summary(Min3vio_only)



majvio_only <- lmer(socvio_only~ MINPCTX + YouthRatio.y+ logland+ ln_population + ln_GDP + popvote_im + church + relimp + majority + bur_qua +  time + lag_MMXXd  +  lag_MMXXmean_int + lag_lxxd_int +lag_lxxmean_int +(1 | country), final_data)

summary (majvio_only)


# 9 - Diagnostics  (Figure A8 and A9)

# Residual normality

# Minority (Model M3 in the appendix)


resmin <- hlm_resid(Min3,include.ls = FALSE)
res2min <-  hlm_resid(Min3, level = "country", include.ls = FALSE)

car::qqPlot(res2min$.ranef.intercept, ylab = "Random effect residuals") 
car::qqPlot(resmin$.resid, envelope = list(level = 0.99), ylab = "Conditional residuals")
car::qqPlot(resmin$.mar.resid, envelope = list(level = 0.99), ylab = "Marginal residuals")


# Majority (Model M7)

res <- hlm_resid(Maj3,include.ls = FALSE)
res2 <-  hlm_resid(Maj3, level = "country", include.ls = FALSE)

car::qqPlot(res2$.ranef.intercept,ylim = c(-5,5),col.lines = "orange",ylab = "Random effect residuals") 
car::qqPlot(res$.resid, envelope = list(level = 0.99),col.lines = "orange",  ylab = "Conditional residuals") 
car::qqPlot(res$.mar.resid, envelope = list(level = 0.99), col.lines = "orange",  ylab = "Marginal residuals") 

# Multicollineartiy (Table A5. Vif Scores)
 
# Minortiy

library(car)

vif(Min3)

vif(Maj3)

# 10 - Bivariate relationships (Table A7)

bivariate_1 <- lmer(WZMIN2MAJX ~  lag_MMXXd  +  lag_MMXXmean_int + time +  (1 | country), final_data)
bivariate_2 <- lmer(WZMIN2MAJX ~  time + lag_lxxd_int  +  lag_lxxmean_int +  (1 | country), final_data)

bivariate_3 <- lmer(soc_vioc_int ~  lag_MMXXd  +  lag_MMXXmean_int + time +  (1 | country), final_data)
bivariate_4 <- lmer(soc_vioc_int ~   time + lag_lxxd_int  +  lag_lxxmean_int  +  (1 | country), final_data)

summary(bivariate_4 )


# - Standardized coefficient plots (Figures A6 and A7)

library(sjPlot)


plot_model(Min3, type = "std") #(Figure A6)

plot_model (Maj3, type = "std") #(Figure A7)


```



*Individual item Analyses (Table A11)*

```{r}
# Muslims 


summary(lmer(WZMIN2MAJX~  lag_MMX06X +MINPCTX+ gini_disp +YouthRatio.y + logland +ln_population + ln_GDP + time  +lag_lxxd_int + lag_lxxmean_int+ (1 | country), final_data)) # 0.351 t value :2.307


summary(lmer(WZMIN2MAJX~  lag_MMX16X +MINPCTX+ gini_disp +YouthRatio.y + logland +ln_population + ln_GDP + time  +lag_lxxd_int + lag_lxxmean_int+ (1 | country), final_data)) #  5.405e-01 # 11.334
 

summary(lmer(WZMIN2MAJX~  lag_MMX36X +MINPCTX+ gini_disp +YouthRatio.y + logland +ln_population + ln_GDP + time  +lag_lxxd_int + lag_lxxmean_int+ (1 | country), final_data)) # 0.398 # 3.924

#Majority

summary(lmer(soc_vioc_int~ lag_MMX01X+  MINPCTX + YouthRatio.y+ logland+ ln_population + ln_GDP + popvote_im + church + relimp + majority + bur_qua + lag_lxxd_int + lag_lxxmean_int+ time + (1 | country), final_data)) #   1.522  #2.534

summary(lmer(soc_vioc_int~  lag_MMX05X+  MINPCTX + YouthRatio.y+ logland+ ln_population + ln_GDP + popvote_im + church + relimp + majority + bur_qua + lag_lxxd_int + lag_lxxmean_int+ time + (1 | country), final_data)) #   3.009  # 7.99

summary(lmer(soc_vioc_int~  lag_MMX10X+  MINPCTX + YouthRatio.y+ logland+ ln_population + ln_GDP + popvote_im + church + relimp + majority + bur_qua + lag_lxxd_int + lag_lxxmean_int+ time + (1 | country), final_data)) #7.020 #3.450

summary(lmer(soc_vioc_int~  lag_MMX16X+  MINPCTX + YouthRatio.y+ logland+ ln_population + ln_GDP + popvote_im + church + relimp + majority + bur_qua + lag_lxxd_int + lag_lxxmean_int+ time + (1 | country), final_data)) #  0.7465  #3.714

summary(lmer(soc_vioc_int~  lag_MMX31X+  MINPCTX + YouthRatio.y+ logland+ ln_population + ln_GDP + popvote_im + church + relimp + majority + bur_qua + lag_lxxd_int + lag_lxxmean_int+ time + (1 | country), final_data)) #   2.9300 #3.528

summary(lmer(soc_vioc_int~  lag_MMX32X+  MINPCTX + YouthRatio.y+ logland+ ln_population + ln_GDP + popvote_im + church + relimp + majority + bur_qua + lag_lxxd_int + lag_lxxmean_int+ time + (1 | country), final_data)) #  0.7153 # 4.186

summary(lmer(soc_vioc_int~  lag_MMX35X+  MINPCTX + YouthRatio.y+ logland+ ln_population + ln_GDP + popvote_im + church + relimp + majority + bur_qua + lag_lxxd_int + lag_lxxmean_int+ time + (1 | country), final_data)) # 3.5103 # 3.450


```


*Augmented Synthetic control(ASCM)* 


#Majority Violence

```{r}

library(augsynth)
final_data <- final_data %>% ungroup()


## Exclude countries that potentially have the same treatment

condition1 <- !(final_data$country %in% c("France", "Denmark", "United States of America"))
final_data_aug1 <- final_data[condition1, ]


# ASCM without covariates ( Table A10 in the appendix)

maj_aug_0 <- augsynth(soc_vioc_int ~ D ,iso3, year, final_data_aug1 , progfunc = "ridge", scm = T)
summary(maj_aug_0, inf_type = "jackknife")

# ASCM with covariates (Table 4- A1 model in the paper)

maj_aug <- augsynth(soc_vioc_int ~ D | majority + ln_population +ln_GDP + church + relimp+MINPCTX + popvote_im + YouthRatio.y +MMXX + logland + lag_socvioc ,iso3, year, final_data_aug1 , progfunc = "ridge", scm = T)

# Summary Results

summary(maj_aug, inf = T, inf_type = "jackknife")

# Model Fit
maj_aug$scaled_l2_imbalance
maj_aug$weights
maj_aug$scaled_covariate_l2_imbalance


# Graph Weigths (Figure A12- left)

country_codes <-  c("AUS", "AUT", "BEL", "CAN", "CYP", "FIN", "DEU","GRC", "ISL", "IRL","ITA", "LIE", "LUX","MLT", "NLD","NZL","NOR","PRT","ESP",
"SWE","CHE")

majx_values <- as.data.frame(maj_aug$weights)

majweight_data <- cbind(majx_values, country_codes)


ggplot(majweight_data, aes(x = V1, y = factor(country_codes))) +
geom_point(color = "blue", size = 3) +
     labs(x = "Weights", y = "", title = "") +
     theme_minimal()

# Covariate Balance table (Table A8/model A1)

cov_synth <- maj_aug$data$Z[c(1:21), ] 

maj_aug$data$Z

# Multiply with weights 
cov_synth <- cov_synth * maj_aug$weights [,1]

# sum  all columns 

cov_synth <- colSums(cov_synth)
cov_synth <- as.data.frame(cov_synth)
cov_UK <- maj_aug$data$Z[22, ]

cov_synth_dat <- cbind(cov_synth, cov_UK)


# Gap plot (Figure 4-a)


 plot(maj_aug,inf_type = "jackknife") + labs(x = "Year", y = "Difference in Majority Violence" ) + geom_vline(xintercept = 2006, linetype = "dashed", color = "red") +
        annotate("segment", x = 2006, xend = 2008, y = 4, yend = 4,
           colour = "blue", size = 0.5, arrow = arrow()) + 
   annotate("text", x = 2010, y = 4, label = "made public", colour = "blue", vjust = 0)
 
# Path plot (Figure 4-b)
y0plot1 <- maj_aug$data$synth_data$Y0plot %*% maj_aug$weights

#Ylim 
  Y.max <-  max(c(y0plot1,maj_aug$data$synth_data$Y1plot))
  Y.min <-  min(c(y0plot1,maj_aug$data$synth_data$Y1plot))
  Ylim <- c((Y.min - .3*Y.min ), (.3*Y.max + Y.max))
 
plot(
                     maj_aug$data$time
                     ,maj_aug$data$synth_data$Y1plot,
                     t="l",
                     col="black",
                     lwd=2, 
                     ylim = Ylim,
                     ylab = "Violence by Majority Groups",
                     xlab ="Year")
                     
                     lines(
                            maj_aug$data$time,
                           y0plot1 ,
                           col="black",
                           lty="dashed",
                           lwd=2,
                           cex=4/5,
                           type = 'l'
                           )

legend("topright", legend = c("Synthetic UK",  "UK"), lty = c(2,1),  pch = 16, col = "black")

abline(v = c(2003, 2006), col = c("black", "red"), lty = 2)
                     
```


*Muslim violence*

```{r}
# without covariates  (Table A10)

min_aug_0 <- augsynth(WZMIN2MAJX ~ D, country, year, final_data_aug1, progfunc = "ridge", scm = T)
summary(min_aug_0, inf_type ="jackknife")

# With covariates (Table 4 model A1)
min_aug <- augsynth(WZMIN2MAJX ~ D | ln_population + ln_GDP + MINPCTX + YouthRatio.y +gini_disp + logland+ MMXX + lag_WZMIN2MAJX ,
country, year, final_data_aug1, progfunc = "ridge", scm = T)

# Summary results

summary(min_aug, inf = T, inf_type = "jackknife")

# Model fit
min_aug$scaled_l2_imbalance
min_aug$weights
min_aug$scaled_covariate_l2_imbalance

# create weight graph (Figure A13)

country_codes <-  c("AUS", "AUT", "BEL", "CAN", "CYP", "FIN", "DEU","GRC", "ISL", "IRL","ITA", "LIE", "LUX","MLT", "NLD","NZL","NOR","PRT","ESP",
"SWE","CHE")

minx_values <- as.data.frame(min_aug$weights)

minweight_data <- cbind(minx_values, country_codes)


ggplot(minweight_data, aes(x = V1, y = factor(country_codes))) +
geom_point(color = "blue", size = 3) +
     labs(x = "Weights", y = "", title = "") +
     theme_minimal()


# Covariate Balance Table (Table A9/model A1)

cov_synth_min <- min_aug$data$Z[c(1:21), ] 

# Multiply with weights 

cov_synth_min <- cov_synth_min * min_aug$weights [,1]

# sum all columns 

cov_synth_min <- colSums(cov_synth_min)
cov_synth_min <- as.data.frame(cov_synth_min)
cov_UK_min <- min_aug$data$Z[22, ]

# Compare weights of synthetic Uk and UK

cov_synth_dat_min <- cbind(cov_synth_min, cov_UK_min)


# Gap plot (Figure 4-c)

plot(min_aug,  inf_type = "jackknife") + labs(x = "Year", y = "Difference in Minority Violence") +geom_vline(xintercept = 2006, linetype = "dashed", color = "red") +         annotate("segment", x = 2006, xend = 2008, y = 1.5, yend = 1.5,
           colour = "blue", size = 0.5, arrow = arrow()) + 
   annotate("text", x = 2010, y = 1.5, label = "made public", colour = "blue", vjust = 0)


# Path Plot (Figure 4-d)

y0plot1_min <- min_aug$data$synth_data$Y0plot %*% min_aug$weights

  Y2.max <-  max(c(y0plot1_min,min_aug$data$synth_data$Y1plot))
  Y2.min <-  min(c(y0plot1_min,min_aug$data$synth_data$Y1plot))
  Ylim2 <- c((Y2.min - .3*Y2.min ), (.3*Y2.max + Y2.max))
 
  plot(
                      min_aug$data$time
                     ,min_aug$data$synth_data$Y1plot,
                     t="l",
                     col="black",
                     lwd=2, 
                     ylim = Ylim2,
                     xlab= "Year",
                     ylab = "Muslim Minority Violence")
                     
                     lines(
                            min_aug$data$time,
                           y0plot1_min ,
                           col="black",
                           lty="dashed",
                           lwd=2,
                           cex=4/5
                           )


legend("topright", legend = c("Synthetic UK",  "UK"), lty = c(2,1),  pch = 16, col = "black")

abline(v = c(2003, 2006), col = c("black", "red"), lty = 2)
                     


```

*In time Placebos*

```{r}
### Majority violence 

# 1999 (Table A10)

placebo_df <- final_data_aug1
placebo_df <- filter(placebo_df,  placebo_df$year < 2004)
placebo_df$D3 <- ifelse(placebo_df$country == "United Kingdom" & placebo_df$year >= 1999, 1, 0)


maj_aug_p1 <- augsynth(soc_vioc_int ~ D3 | majority + ln_population + ln_GDP + church + relimp +MINPCTX + popvote_im + bur_qua  + popvote_im+YouthRatio.y+MMXX + lag_socvioc + logland,iso3, year, placebo_df, progfunc = "ridge", scm = T)

summary(maj_aug_p1, inf = T, inf_type = "jackknife")

# Gap Plot (Figure A15)

plot(maj_aug_p1, inf_type = "jackknife") + labs(x = "Year", y = "Difference in Majority Violence")


# 1996 (Table A10)

placebo_df$D4 <- ifelse(placebo_df$country == "United Kingdom" & placebo_df$year >= 1996, 1, 0)


maj_aug_p2 <-  augsynth(soc_vioc_int ~ D4 | majority + ln_population + ln_GDP + church + relimp +MINPCTX + popvote_im + bur_qua  + popvote_im+YouthRatio.y +MMXX + lag_socvioc + logland, iso3, year, placebo_df, progfunc = "ridge", scm = T)

summary(maj_aug_p2, inf = T, inf_type = "jackknife")

# Figure A15
plot(maj_aug_p2, inf_type = "jackknife") + labs(x = "Year", y = "Difference in Majority Violence")


### Muslim violence 

# 1999 (Table A10)
min_aug_p2 <-  augsynth(WZMIN2MAJX ~ D3 | ln_population + ln_GDP + MINPCTX + YouthRatio.y +gini_disp + logland+ MMXX + lag_WZMIN2MAJX,
country, year, placebo_df, progfunc = "ridge", scm = T)

summary(min_aug_p2, inf = T, inf_type = "jackknife")
plot(min_aug_p2, inf_type = "jackknife")


#1996 

min_aug_p1 <-  augsynth(WZMIN2MAJX ~ D4| ln_population + ln_GDP + MINPCTX + YouthRatio.y +gini_disp + logland+ MMXX +lag_WZMIN2MAJX,
country, year, placebo_df, progfunc = "ridge", scm = T)

summary(min_aug_p1, inf = T, inf_type = "jackknife")
plot(min_aug_p2, inf_type = "jackknife")



```

*Exclude  additional control countries for robustness (Model A2)*


```{r}
### Majority violence 

condition2 <- !(final_data$country %in% c("France", "Denmark", "United States of America", "Australia", "New Zealand", "Switzerland", "Netherlands"))

filtered_final_data <- final_data[condition2, ]


maj_augfilt <- augsynth(soc_vioc_int ~ D | majority + ln_population + ln_GDP + church + relimp+MINPCTX + bur_qua  + popvote_im + YouthRatio.y+MMXX + logland + lag_socvioc, iso3, year, filtered_final_data, progfunc = "ridge", scm = T)

# Table 4-A2
summary(maj_augfilt, inf_type ="jackknife") 

# Figure A14-a

plot(maj_augfilt,  inf_type = "jackknife")

# Model Fit
maj_augfilt$scaled_l2_imbalance
maj_augfilt$weights

# Graph Weights (Figure A12- right)

country_codes_filt <-  c( "AUT", "BEL", "CAN", "CYP","FIN", "DEU","GRC", "ISL", "IRL","ITA", "LIE", "LUX","MLT", "NOR","PRT","ESP",
"SWE")

majx_values_filt <- as.data.frame(maj_augfilt$weights)

majweight_data_filt <- cbind(majx_values_filt, country_codes_filt)


ggplot(majweight_data_filt, aes(x = V1, y = factor(country_codes_filt))) +
geom_point(color = "orange", size = 3) +
     labs(x = "Weights", y = "", title = "") +
     theme_minimal()


# Covariate Balance Table (table A8)

cov_synth_filt <- maj_augfilt$data$Z[c(1:17), ] 

# Multiply with weights 
cov_synth_filt <- cov_synth_filt * maj_augfilt$weights [,1]

# sum all columns 

cov_synth_filt <- colSums(cov_synth_filt)
cov_synth_filt <- as.data.frame(cov_synth_filt)
cov_UK_filt <- maj_augfilt$data$Z[18, ]

# Combine and compare covariates

cov_synth_dat_filt <- cbind(cov_synth_filt, cov_UK_filt)

# Gap Plot - Figure A14-a


plot(maj_augfilt,  inf_type = "jackknife") +  labs(x = "Year", y = "Difference in Majority Violence")


###  Minority Violence 

min_augfilt <- augsynth(WZMIN2MAJX ~ D | ln_population + ln_GDP + MINPCTX  + YouthRatio.y +gini_disp + logland+ MMXX +lag_WZMIN2MAJX,
iso3, year, filtered_final_data, progfunc = "ridge", scm = T)

# Table 4- model A2

summary(min_augfilt, inf = T, inf_type = "jackknife")

# Model Fit
min_augfilt$scaled_l2_imbalance
min_augfilt$weights

# Graph Weights  (Figure A13-right)

country_codes_filt <-  c("AUT", "BEL", "CAN", "CYP","FIN", "DEU","GRC", "ISL", "IRL","ITA", "LIE", "LUX","MLT", "NOR","PRT","ESP","SWE")

minx_values_filt <- as.data.frame(min_augfilt$weights)

minweight_data_filt <- cbind(minx_values_filt, country_codes_filt)


ggplot(minweight_data_filt, aes(x = V1, y = factor(country_codes_filt))) +
geom_point(color = "orange", size = 3) +
     labs(x = "Weights", y = "", title = "") +
     theme_minimal()


##  Covariate Balance Table (Table A9)

cov_synth_filt_min <- min_augfilt$data$Z[c(1:17), ] 

# Multiply with weights 
cov_synth_filt_min <- cov_synth_filt_min * min_augfilt$weights [,1]

# sum all columns 

cov_synth_filt_min <- colSums(cov_synth_filt_min)
cov_synth_filt_min <- as.data.frame(cov_synth_filt_min)
cov_UK_filt_min <- min_augfilt$data$Z[18, ]

#  Combine and compare covariates

cov_synth_dat_filt_min <- cbind(cov_synth_filt_min, cov_UK_filt_min)

#  Gap Plot  (Figure A14-b)

plot(min_augfilt,  inf_type = "jackknife") + labs(x = "Year", y = "Difference in Minority Violence")

# Model including transnational terrorism (as footnote 16 in the paper)

# ASCM with covariates 

maj_augx <- augsynth(soc_vioc_int ~ D | majority + ln_population +ln_GDP + church + relimp+MINPCTX + popvote_im + YouthRatio.y +MMXX + logland + lag_socvioc+ lag(transnationalterrorism) ,iso3 , year, final_data_aug1 , progfunc = "ridge", scm = T)

summary(maj_augx, inf = T, inf_type = "jackknife")

```


